# -*- coding: utf-8 -*-
"""
@author:xuyuntao
@time:2021/1/14:8:07
@email:xuyuntao@189.cn
"""
import numpy
import numpy.fft as nfft
import scipy.signal as scisignal
import csv
from Tools import runTimeCounter
import time

showTime=False

class AnalogModulationDetect(object):
    def __init__(self,signal=None,sample_freq=1e3,\
                 t_a=1,t_gamma_max=500,t_sigma_dp=0.5,t_P=0.6,\
                 Amp_pass=1,Amp_stop=30):
        if (type(signal)==numpy.ndarray and len(signal.shape)==1):
            self.signal=signal.astype(numpy.float64)
        else:
            if signal==None:
                self.signal=numpy.zeros(1000,dtype=numpy.float64)
            else:
                raise ValueError
        self.f_sample=sample_freq # 采样频率
        self.T_sample = 1 / self.f_sample # 采样周期
        self.t_a=t_a
        self.t_gamma_max=t_gamma_max # 5.8  # gamma_max的特征门限值
        self.t_sigma_dp=t_sigma_dp # sigma_dp的特征门限值
        self.t_P=t_P # P的特征门限值
        self.time_series=numpy.linspace(0, self.signal.size - 1, self.signal.size) * self.T_sample
        self.Amp_pass=Amp_pass # Ap：通带最大衰减
        self.Amp_stop=Amp_stop # As：阻带最小衰减
        self.modulationName=["AM", "DSBSC", "SSB", "FM"]
        self.modulationType=-1
        self._signalFFT()

    @runTimeCounter(showTime)
    def _signalFFT(self):
        """计算信号频谱"""
        self.signal_fft = nfft.fft(self.signal)
        return self.signal_fft

    @runTimeCounter(showTime)
    def _getSignalFreq(self):
        """信号载频估计"""
        # 基于过零点检测估计频率
        signal_augment = numpy.zeros(self.signal.size + 1)
        signal_augment[1:self.signal.size + 1] = self.signal
        signal_forword = signal_augment.copy()
        signal_augment[0:self.signal.size] = self.signal
        signal_augment[-1] = 0
        signal_after = signal_augment.copy()
        zero_crossing = (signal_forword * signal_after) < 0
        zero_crossing_num = numpy.sum(zero_crossing[2:self.signal.size])

        last_zero_2 = numpy.where(zero_crossing[0:self.signal.size])[0][-1]
        last_zero_2_value = self.signal[last_zero_2]
        last_zero_1 = last_zero_2 - 1
        last_zero_1_value = self.signal[last_zero_1]

        first_zero_2 = numpy.where(zero_crossing[0:self.signal.size])[0][0]
        first_zero_2_value = self.signal[first_zero_2]
        first_zero_1 = first_zero_2 - 1
        first_zero_1_value = self.signal[first_zero_1]

        # 计算载频
        z_last = last_zero_2 - last_zero_2_value * (last_zero_1_value - last_zero_2_value)  # 最后一个过零点区域斜率
        z_first = first_zero_2 - first_zero_2_value * (first_zero_1_value - first_zero_2_value)

        self.f_carrier = 0.5 * self.f_sample * (zero_crossing_num - 1) / (z_last - z_first)
        self.f_Res = self.f_sample / self.signal.size  # 频率分辨率
        self.f_carrier_fft_index = round(self.f_carrier / self.f_Res)  # 中心频率f0在X(k)中的位置
        self.kc=self.f_carrier_fft_index

        return self.f_carrier

    @runTimeCounter(showTime)
    def _designFilter(self):
        """设计滤波器传输系数"""
        f_pass = round(self.f_carrier)  # fp：滤波器的通带频率
        f_stop = round(0.9 * self.f_carrier)  # fst: 滤波器的截止频率
        f_pass_digital = 2 * f_pass / self.f_sample
        f_stop_digital = 2 * f_stop / self.f_sample
        self._filterSteps, self._omega_n = scisignal.buttord(f_pass_digital, f_stop_digital,\
                                                 self.Amp_pass, self.Amp_stop)
        self.filter_b, self.filter_a = scisignal.butter(self._filterSteps, self._omega_n,\
                                                        "low", output="ba")

    @runTimeCounter(showTime)
    def _quadMixLowPass(self):
        """正交解调+低通滤波"""
        # 混频
        w_signal = 2 * numpy.pi * self.f_carrier  # 模拟角频率
        signal_inPhase = self.signal * 2 * numpy.cos(w_signal * self.time_series)
        signal_quadPhase = self.signal * -2 * numpy.sin(w_signal * self.time_series)

        # 低通滤波
        try:
            self.filter_a
            self.filter_b
            # print("Filter A = ",self.filter_a)
            # print("Filter B = ",self.filter_b)
        except NameError:
            self._designFilter()
            # print("Filter A = ",self.filter_a)
            # print("Filter B = ",self.filter_b)
        finally:
            self.signal_inPhase_lowPass = scisignal.filtfilt(self.filter_b, self.filter_a, signal_inPhase)
            self.signal_quadPhase_lowPass = scisignal.filtfilt(self.filter_b, self.filter_a, signal_quadPhase)

        return self.signal_inPhase_lowPass,self.signal_quadPhase_lowPass

    @runTimeCounter(showTime)
    def _getP(self):
        """计算P：谱对称性"""
        P_L = numpy.sum(numpy.abs(self.signal_fft[0:self.kc]) ** 2)
        P_U = numpy.sum(numpy.abs(self.signal_fft[self.kc + 1:2 * self.kc + 1]) ** 2)
        self.P = (P_L - P_U) / (P_L + P_U)
        return self.P

    @runTimeCounter(showTime)
    def _getCharacters(self):
        """计算gamma_max和sigma_dp"""
        # ---------计算各个特征值---------
        phi = numpy.arctan2(self.signal_inPhase_lowPass, self.signal_quadPhase_lowPass)
        phi_avg = numpy.mean(phi)  # 平均值
        phi_NL = phi - phi_avg
        amp = numpy.sqrt(self.signal_inPhase_lowPass ** 2 + self.signal_quadPhase_lowPass ** 2)  # 幅度
        amp_avg = numpy.mean(amp)  # 平均值
        amp_n = amp / amp_avg  # 用平均值对瞬时幅度进行归一化
        amp_cn = amp_n - 1

        # ---------计算特征值gamma_max：归一化零中心瞬时幅度之谱密度的最大值---------
        self.gamma_max = numpy.max(numpy.abs(nfft.fft(amp_cn))) ** 2 / self.signal.size

        # ---------计算特征值sigma_dp：非弱信号段零中心瞬时相位非线性分量的方差---------
        judge_A_n = amp_n > self.t_a
        self.sigma_dp = numpy.sqrt(numpy.sum(phi_NL[judge_A_n] ** 2) / numpy.sum(judge_A_n) -\
                                   (numpy.sum(phi_NL[judge_A_n]) / numpy.sum(judge_A_n)) ** 2)

        return self.gamma_max,self.sigma_dp

    @runTimeCounter(showTime)
    def _detectModulation(self):
        """识别调制方式"""
        if (self.sigma_dp < self.t_sigma_dp):
            self.modulationType = 0  # 常规调幅AM
        else:
            if (abs(self.P) >= self.t_P):
                self.modulationType = 2  # 为SSB, 本实验中，其为下边带LSB
            else:
                if (self.gamma_max < self.t_gamma_max):
                    self.modulationType = 3  # 为FM
                else:
                    self.modulationType = 1  # 为DSBSC

    @runTimeCounter(showTime)
    def _getSpectrumWidth(self):
        """计算信号功率、90%功率和、信号带宽"""
        if (self.modulationType == 0):  # 若为常规调幅AM
            # 对于AM调制方式，由于其在载波频率处有一个大的冲激，其占据信号功率的绝大部分，
            # 因此，对于载频点在FFT中的位置fcn，必须计算完全准确，否则将给功率计算带来巨大误差，
            # 因此，为避免f0的少许误差使得fcn的值发生偏差，在上面已识别出是AM调制方式的前提下，
            # 下面重新计算fcn的位置
            self.f_carrier_fft_index = numpy.where(numpy.abs(self.signal_fft) ==\
                                                   numpy.max(numpy.abs(self.signal_fft)))[0][0]
            # print("AM位置为index",self.f_carrier_fft_index)
        # ---------计算总的功率---------
        self.TotalPower = numpy.sum(numpy.abs(self.signal_fft[0:round(self.signal.size / 2)]) ** 2)  # 总的功率
        if (self.f_carrier_fft_index) < round(self.signal.size / 2):
            self.TotalPower -= numpy.abs(self.signal_fft[self.f_carrier_fft_index]) ** 2

        # ---------计算90%功率之和psum---------
        startTime=time.time()
        temp_value = numpy.abs(self.signal_fft[self.f_carrier_fft_index])  # 载频处的值在kc位置
        self.psum = 0
        if (self.modulationType == 0):
            self.psum = 0  # 若是AM调制方式，psum不应包括载频处的功率，故psum初值为0
        else:
            self.psum = temp_value ** 2  # 初值为载频处的功率
        psum_stop = 0
        # for i=kc:-1:1 # 从载频处开始计算功率之和
        for _ in range(self.f_carrier_fft_index-1,0,-1):
            if (self.modulationType == 2):  # 若是单边带(下边带)调制，频谱不对称
                self.psum += numpy.abs(self.signal_fft[_]) ** 2
            else:  # 若是其他调制，频谱对称
                self.psum += numpy.abs(self.signal_fft[_]) ** 2  # 载频左边的功率
                self.psum += numpy.abs(self.signal_fft[self.f_carrier_fft_index + (self.f_carrier_fft_index - _)]) ** 2  # 载频右边的功率
            if (self.psum / self.TotalPower >= 0.9):  # 信号带宽按总功率的90%来确定
                psum_stop = _
                break
        endTime=time.time()
        # print("----计算90%功率之和psum使用时间为{0:.4f}秒----".format(endTime-startTime))

        # ---------计算频谱宽度SpectrumWidth---------
        if (self.modulationType == 2):
            self.spectrumWidth = (self.f_carrier_fft_index - psum_stop) * self.f_Res
        else:
            self.spectrumWidth = (self.f_carrier_fft_index - psum_stop) * 2 * self.f_Res

        return self.TotalPower,self.psum,self.spectrumWidth

    def updateSignal(self,signal):
        """更改signal"""
        self.signal=signal
        self.time_series = numpy.linspace(0, self.signal.size - 1, self.signal.size) * self.T_sample
        self._signalFFT()
        return self.signal

    def getAllOnce(self):
        """一次算完"""
        startTime=time.time()
        self._getSignalFreq()
        self._designFilter()
        self._quadMixLowPass()
        self._getP()
        self._getCharacters()
        self._detectModulation()
        self._getSpectrumWidth()
        self.printResult()
        endTime=time.time()
        print("----总计算时间为{0:.4f}秒----\n".format(endTime-startTime))

    def printResult(self):
        """显示结果"""
        print("P = {3}\ngamma_max = {4}\nsigma_dp = {5}\n调制方式：{0}\n带宽：{1:.4f}kHz\n载波频率：{2}Hz".\
              format(self.modulationName[self.modulationType],\
                     self.spectrumWidth/1000, self.f_carrier, self.P,\
                     self.gamma_max, self.sigma_dp))

if __name__=="__main__":
    dataPath="AnalogModulationSignal"
    signal=numpy.zeros(10,dtype=numpy.float32)
    DAM = AnalogModulationDetect(signal, 8e7)
    for fileNum in range(4):
        startTime=time.time()
        # signal=numpy.loadtxt(dataPath+"\\"+"data{}.csv".format(fileNum+1), dtype=numpy.float32, delimiter=',')
        data=[]
        with open(dataPath+"\\"+"data{}.csv".format(fileNum+1),"r") as f:
            csvRead=csv.reader(f)
            for _ in csvRead:
                data.append(_)
        signal=numpy.array(data,dtype=numpy.float32).reshape(-1)
        print("--------------------文件：{}--------------------".format("data{}.csv".format(fileNum + 1)))
        DAM.updateSignal(signal)
        DAM._getSignalFreq()
        DAM._designFilter()
        DAM._quadMixLowPass()
        DAM._getP()
        DAM._getCharacters()
        DAM._detectModulation()
        DAM._getSpectrumWidth()
        DAM.printResult()
        endTime=time.time()
        print("----总计算时间为{0:.4f}秒----\n".format(endTime-startTime))


